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The value measured in the amorphous structure with the same chemical composition is often con¬ 
sidered as a lower bound for the thermal conductivity of any material: the heat carriers are strongly 
scattered by disorder, and their lifetimes reach the minimum time scale of thermal vibrations. An 
appropriate design at the nano-scale, however, may allow one to reduce the thermal conductivity 
even below the amorphous limit. In the present contribution, using molecular-dynamics simula¬ 
tion and the Green-Kubo formulation, we study systematically the thermal conductivity of layered 
phononic materials (superlattices), by tuning different parameters that can characterize such struc¬ 
tures. We discover that the key to reach a lower-than-amorphous thermal conductivity is to block 
almost completely the propagation of the heat carriers, the superlattice phonons. We demonstrate 
that a large mass difference in the two intercalated layers, or weakened interactions across the in¬ 
terface between layers result in materials with very low thermal conductivity, below the values of 
the corresponding amorphous counterparts. 

I. INTRODUCTION 

Materials with low thermal conductivity, k, are em¬ 
ployed in many modern technologies, such as thermal 
management in electronic devices or thermoelectric en¬ 
ergy conversion [1-3]. In general, low values of k are ob¬ 
served in disordered solids [4], including topologically dis¬ 
ordered systems (glasses) and crystalline solids with size 
or mass disorder (disordered alloys) [5-9]. This behaviour 
can be rationalized by considering the phenomenological 
kinetic theory expression [10] n = (1/3)C'u^t, which re¬ 
lates average velocity, v, and lifetime, r, (and therefore 
the mean free path £ = vt) of phonons to n {C is the spe¬ 
cific heat per unit volume). In good crystals, phonons 
lifetime is primarily controlled by anharmonic interac¬ 
tions. In contrast, in disordered solids, the disorder (or 
the elastic heterogeneity [11]) reduces r (or £) and, as a 
result, K. 

In early experimental investigations [5, 6], Cahill 
et al. have studied the disordered alloys, e.g., 

(KBr)i_a,(KCN)a;, and shown that k can be reduced to 
the glass value by controlling the relative composition 
X. In our works [8, 9] we in turn demonstrated that, 
in size-disordered crystal, k progressively decreases with 
increasing size mismatch, eventually converging to the 
corresponding glass value. When this limit is reached, r 
is comparable to the time scale of thermal vibrations (£ 
to the particle size), i.e., to the minimum time (length) 
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scale [8, 9]. Heat propagation can therefore be described 
as a random walk of vibrational energies [5, 6], or in 
terms of non-propagating delocalized modes, the diffu- 
sons [7]. For this reason, the value in the glass is gener¬ 
ally considered as a lower bound for k of materials with 
homogeneous chemical composition [5, 6]. 

A crucial issue [4] is whether thermal conductivity 
can be lowered below the glass limit through nanoscale 
phononic design [3, 12]. This possibility would allow to 
devise (meta-)materials which are excellent thermal in¬ 
sulators while preserving good electronic properties, as 
needed in many applications [1-3]. The most popular 
design to reach this goal is that of a lamellar super¬ 
lattice [13-17], often composed of two chemically dif¬ 
ferent intercalated layers, e.g., Si-Ge [13, 14] or GaAs- 
AlAs [15, 16] (see also Fig. 1). In a superlattice, the 
thermal conductivity tensor is anisotropic, with the cross¬ 
plane component, kcP: usually lower than the in-plane 
value, Kip [18, 19]. In recent experiments [20-22], ultra- 
low values of kcp i suggested to be smaller than the amor¬ 
phous limit, were measured. In particular, Gostescu et 
al. [20] demonstrated that the presence of a high-density 
of interfaces decreases kqp of W-AI2O3 nanolaminates, 
below that of the amorphous AI2O3. An experiment by 
Ghiritescu et al. [21] achieved ultra-low thermal conduc¬ 
tivity in layered WSe2 crystals, by disordering the crys¬ 
talline WSe2 sheets. Finally, Pernot et al. [22] also ob¬ 
served very low values of kcPj below that of amorphous 
Si, in Ge nanodots multi-layers separated by Si crystals. 

Although the above works have demonstrated very low 
values of n in superlattice systems, we note that these 
have not been systematically compared to the values as¬ 
sumed in the glasses with exactly the same chemical com¬ 
position. Also, a general framework to rationalize in a 
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TABLE I. The investigated superlattice structures. Details of the three superlattice systems investigated in this work. 
They are based on the FCC-crystal lattice structure and are composed of: (SI) two intercalated crystalline layers {A and B) 
formed by sphere particles with different masses iriA and ms; {S2) ordered crystalline layers intercalated to mass-disordered 
alloy layers; and {S3) identical crystalline layers with modified (weakened compared to those intra-layers) interactions across 
the interfaces. The control parameters are the mass ratio ms/rnA in SI, the mass ratio rriBiIrnBi of the disordered alloy layer 
in S2, and the energy scale tAB of the interactions across the interfaces in S3. Number density and temperature were fixed to 
the values p = 1.015 (corresponding to a lattice constant a = 1.58) and T = 10~^, respectively. The quantities presented in 
the table are defined in the main text. In the last column we refer to the figure containing the data relative to the indicated 
system. Additional details about the investigated superlattices and parameters used are given in the Methods section. 


System 

Control Parameter 


ka 

Kb 

KcP 

Bip 

R 

£k 

^glass 

^disorder 

Fig. 

{SI) Mass difference 

niB/mA = 

2 

488.6 

335.4 

397.8 

412.0 

0.5 

398 

5.7 

20.4 

Fig. 2(a) 



4 

625.9 

306.8 

411.8 

466.3 

1.9 

1564 

4.2 

9.9 

Fig. 2(b) 



8 

843.8 

291.1 

432.8 

567.5 

- 

- 

3.3 

7.7 

Fig. 2(c) 

{S2) Order-disorder 

mB^lniBi = 

2 

381.6 

20.4 

38.7 

201.0 

- 

- 

5.7 

33.2 

Fig. 5(a) 



4 

381.6 

9.9 

19.3 

195.8 

- 

- 

4.5 

14.3 

Fig. 5(b) 



8 

381.6 

7.7 

15.1 

194.6 

- 

- 

4.0 

8.2 

Fig. 5(c) 

{S3) Weak interface 

^AB = 

0.5 

587.3 

587.3 

- 

- 

- 

- 

10.6 

- 

Fig. 7(a) 



0.1 

587.3 

587.3 

- 

- 

- 

- 

10.6 

- 

Fig. 7(b) 


coherent single picture all these observations is, to the 
best of our knowledge, still lacking. 

In this work, we address these two issues. Building 
on the comparison of the superlattice with the corre¬ 
sponding amorphous structure, we clarify the mecha¬ 
nisms allowing for ultra-low thermal conductivity in the 
former. We have studied by computer simulation a nu¬ 
merical model that allows one to exactly compare ordered 
and disordered systems with identical chemical composi¬ 
tion and access detailed information on the entire normal 
modes spectrum, providing, as a consequence, a complete 
understanding of the heat transfer process. As the life¬ 
time of heat carriers is already minimum in glasses [8, 9], 
we demonstrate that the key to even lower thermal con¬ 
ductivities is to suppress their propagation across the in¬ 
terfaces between the constituent layers. 

More in details, we have focused on three distinct de¬ 
sign principles for super lattices, mimicking similar con¬ 
figurations actually employed in experiments. These are 
based on the face-centered-cubic (FCC) lattice structure, 
and are composed of: {SI) two intercalated crystalline 
layers formed by sphere particles with different masses; 
{S2) ordered crystalline layers intercalated to mass- 
disordered alloy layers; and {S3) identical crystalline lay¬ 
ers with modified (weakened) interactions across the in¬ 
terfaces (see the Methods section and Table I). We 
show that a large mass difference between layers {SI) 
and weakened interactions between layers {S3) efficiently 
obstruct the propagation of phonons, resulting in a very 
large reduction of the superlattice thermal conductivity, 
even below the values pertaining to the glass phases with 
identical composition. Based on our results, we conclude 
with a discussion of the optimal strategy to follow to¬ 
wards very low thermal conductivity materials. 

In Fig. I we show a schematic illustration of a super¬ 
lattice composed by two intercalated layers, A and B, 


both of thickness w/2. The competition between two 
length scales, the repetition period of the superlattice, 
w, and the mean free path of the superlattice phonons, 
£, determines the coherent or incoherent character of 
phonon transport, as described in [23-25] and demon¬ 
strated by numerical simulations [26-28] and recent ex¬ 
periments [29]. 

For w > £, the incoherent phonon transport is in¬ 
dependent in the different layers, and phonons can 
be effectively treated as particles. In this case, the 
Boltzmann transport equation applies [30, 31], and the 
particle-like phonons are scattered within the layers (in¬ 
ternal resistance) and at the interfaces (interfacial resis¬ 
tance) [32, 33]. The thermal conductivity in the cross¬ 
plane direction can be written as 


KCP = 
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W—¥00 
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Here, ka and kb are the thermal conductivities of ma¬ 
terials A and B, and £k = 2Rk'qp is the Kapitza 
length [34, 35]. R is the interfacial resistance, which ex¬ 
ists even at a perfect interface and depends on the nature 
of the contacting materials (e.g., crystal-crystal, crystal- 
glass) [32, 33]. For w < £k {w > £k), the interfacial 
resistance is relatively large (small) compared to the in¬ 
ternal resistance. Both kcp smd Kip (the in-plane ther¬ 
mal conductivity) increase with w, due to the decrease 
of the interfacial resistance density [30, 31]. In the dif¬ 
fuse limit w —>■ oo, where the interfacial resistance can be 
neglected, kcp and Kip have the upper bounds and 
+ Kb) 12., respectively. 
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When w < £, phonon transport is coherent [23-29], 
and the wave nature of phonons cannot be neglected. 
In this regime, kcp decreases with increasing w, in con¬ 
trast with the incoherent case. The reduction of kqp is 
explained with the emergence of a band gap at the Bril- 
louin zone boundary, due to band-folding [36, 37]: in¬ 
creasing w augments the frequency gap in the dispersion 
relation. This, in turn, decreases the average group ve¬ 
locity V of phonons, finally reducing kcp- Mini-umklapp 
processes [38], occurring at the mini-Brillouin zone, also 
contribute to the reduction of kcp- At the crossover 
length w ^ between the incoherent and the coherent 
transport regimes, kcp assumes a minimum value when 
plotted against w [23-29]. We have encountered this sit¬ 
uation in the case of super lattice SI, as we will see below. 

Details of the structure of the interface between lay¬ 
ers are also known to significantly affect phonon trans¬ 
port [39-49]. It has been reported that interfacial rough¬ 
ness [39-41] or mixing [42, 43] reduce both kcp and kip, 
and can even suppress the coherent nature of phonons, 
with Kcp(ip) increasing monotonously at any w. The in¬ 
terface topology is also an important factor to determine 
the phonon transport [44, 45]. While we will not address 
precisely this situation in detail here, the superlattice S2 
of our study bears some similarities with it. 

Finally, the stiffness of interfacial bondings, which can 
be controlled by applying pressure [46, 47] or tuning 
chemical bonding [48], has significant effects on heat 
transport features, which will be demonstrated by the 
study of the S3 superlattice. 



FIG. 1. Schematic illustration of the considered super¬ 
lattice structures. The investigated superlattice is com¬ 
posed of two FCC-crystalline layers, A (red) and B (green). 
The two layers have identical thickness w/2, where w is the 
replication period. Here, we measure w as the number of 
monolayers of the crystalline lattice, e.g., w = 8 in the dis¬ 
played case. The distance between adjacent monolayers is 
a/2 for theperfect FCC structure we consider, where a is the 
lattice constant. 



W (monolayers) 

FIG. 2. Thermal conductivity in superlattice SI com¬ 
posed of two intercalated crystalline layers with dif¬ 
ferent masses. The cross-plane, kcp, and in-plane, Kip, 
components of thermal conductivity are plotted as functions 
of the repetition period w. The ratio nriB/mA of the masses 
in layers A and B is 2 in panel (a), 4 in (b), and 8 in (c). The 
values k'qb and k)^ of the diffuse limits (w oo), as well 
as those in the glass and the disordered alloy with the same 
constituent species are indicated by the horizontal lines. In 
panels (a) and (b) we also show (dashed black lines), the pre¬ 
diction of Eq. (1) for kcp in the incoherent regime, w > 20, 
with the values of R and Ik included in Table I. The solid 
curve interpolating the kcp data points in the entire w-range 
is a guide for eye. For some values of w, multiple data points 
are shown, calculated by using different system sizes in order 
to exclude the presence of finite system size issues (see the 
Methods section for details on this point). 


II. RESULTS 


In Table I, we present the details of the three su¬ 
perlattice systems studied in this work, with values of 
the important quantities: ka and kb are the thermal 
conductivities of layers A and B, respectively; K^p 
and k[^ are the cross- and in-plane diffuse limits of 
KCP and Kip; R is the interfacial resistance, £k the 
Kapitza length; Kgiass and Kuisorder are the thermal 
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FIG. 3. Vibrational density of states in superlattice SI. Vibrational density of states data for a mass ratio niBlmA = 4, 
with ruA ~ 0.4 and ms = 1.6. In panels (a)-(f) we show the data corresponding to the repetitions period values w = 
2,4,10,20,40,80. For comparison, we also plot gA(B){‘^) for the homogeneous bulk crystal composed by light (heavy) mA(B) 
masses only, together with the data for the glass and the disordered alloy formed by the same constituent species. 


conductivities of the glass and disordered alloy with 
exactly the same composition as the indicated super¬ 
lattice. Thermal conductivities have been estimated by 
molecular-dynamics (MD) simulation and the Green- 
Kubo formulation [50, 51]. The number density and the 
temperature are fixed at p = 1.015 (the corresponding 
crystal lattice constant is a = 1.58) and T = 10“^, 
respectively. Vibrational states were also characterized 
by using a standard normal-modes analysis [52]. Details 
about the systems and the methods used for the sim¬ 
ulation production runs and analysis are given in the 
Methods section. 

SI. Superlattice composed of two intercalated 
crystalline layers with different masses. In Fig. 2 

we show the thermal conductivities, kcp and Kip (sym¬ 
bols), as functions of the replication period, w, for the 
layers mass ratios ms/mA = 2, 4, and 8. The values of 
the diffuse limits Kg’p and as well as those of the glass 
and the disordered alloy constituted by the same species 
(see Table I) are also shown as lines. As expected, the 
relation y/niAKA = y/mBUs holds for the pure materials. 
In the studied rc-range, w = 2 to 40 (monolayers), the 
in-plane value Kip shows a very weak dependence on re, 
as was observed for superlattices with perfect interfaces 
in Refs. [27, 41]. The value of kjp is close to, although 
lower than, k^, indicating that slight in-plane phonon 
scattering at the interface is still active. 

More interestingly, as w increases, the cross-plane 
value Kcp decreases steeply, reaches a minimum value 
at w* ~ 20, and next increases mildly at larger w. This 
w-dependence is consistent with previous predictions [23- 
29], and corresponds to the crossover at w* from coher¬ 
ent to incoherent phonon transport. In the incoherent 
regime, w > 20, from Eq. (1) and the data of kcp (dashed 


line in Fig. 2) we can extract the values of the interfacial 
resistance, R, and the Kapitza length, which are pre¬ 
sented in Table I. Note that for niB/mA = 8 (Fig. 2(c)), 
we do not observe a clear thermal conductivity minimum. 
More precisely, even at the largest value w = 40, kcp is 
still orders of magnitude lower than Kg?, indicating that 
the interfacial resistance R results in a strong reduction 
of KCP in this range of w. Equivalently, the Kapitza 
length £k is signihcantly larger than the maximum pe¬ 
riod w = 40. The data shown in Fig. 2 demonstrate 
that KCP can be indeed lowered below the disordered 
alloy limit for tob/to^ = 2 , and even below the glass 
limit for higher mass heterogeneities, rriB/rriA = 4 and 8. 
These results are consistent with the experimental work 
of Ref. [20], and demonstrate that the interface formed 
between dissimilar materials effectively reduces Kcp. It 
is also worth noting that the thermal conductivity tensor 
is very strongly anisotropic in this case, with kcp <C Kip. 

The vibrational modes of the structure, i.e., the su¬ 
perlattice phonons, are key to understand the above be¬ 
haviour of thermal conductivity. In Fig. 3 we show the vi¬ 
brational density of states (vDOS), g{oj), for mB/mA = 4 
and w = 2 to 80. gA{^) and gB{i^) of the bulk crystals 
of type A and B. The vDOS of the glass and of the 
disordered alloy are also shown for comparison. Note 
that gAj\/"m a = gB^\J"rnB^)/\/‘mB- At small 
w = 2, g{uj) of the superlattice roughly follows that of 
the disordered alloy, implying that the vibrational states 
in the two layers are strongly mixed. In this situation, 
phonons are able to propagate in both the cross- and in¬ 
plane directions. On the other hand, as w increases, gioj) 
generates features increasingly similar to those identify¬ 
ing gA{^) and 5 b(^^); separately. In particular, in the 
low-w region g{uj) follows gB{^) (the heavy crystal B), 
whereas gA{^) (the light crystal A) controls g{io) in the 
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FIG. 4. Vibrational amplitudes of normal modes in superlattice SI. Vibrational amplitudes of eigenvectors, E\ and 
E%, in layers A (light) and B (heavy) for ah normal modes k, plotted as functions of the corresponding eigenfrequency to*. 
E\ and Eg are dehned in Eq. (6). The mass ratio of the two layers is mglmA ~ 4, and the repetitions period values are 
= 2, 4, 10,20,40,80 in panels (a)-(f). The solid line represents the average values {Ea) and (E^) calculated in bins of the 
form u)^ ± 6ui^ 12, with = 0.5. The horizontal dotted lines represent the threshold value Ea = E% = 0.5, the vertical lines 
indicate uj = ~ 13, corresponding to the high frequency edge of ggiuj). 


high-w region. This result indicates that different parts 
of the vibrational spectrum are active in the two lay¬ 
ers, with high(low)-a; modes preferentially excited in the 
light (heavy) layer A (B). In this situation, phonon prop¬ 
agation is largely obstructed in the cross-plane direction, 
leading to the observed large reduction of kcp- We re¬ 
mark that phonons propagate in the in-plane direction 
with few constraints, as shown by the large value of Kip 
close to This implies that phonons, whose propaga¬ 
tions are blocked in the cross-plane direction, are actually 
specularly reflected at the interface and confined in the 
in-plane direction. 

The separation of the vibrational states found in the 
g{uj) becomes more clear when considering the vibra¬ 
tional amplitudes associated with the eigenstates k. In 
Fig. 4 we show the vibrational amplitudes, and Eg 
(Eq. (6)), in the two layers A and B for each mode k, to¬ 
gether with the binned average values (solid lines). Based 
on the relations E^ + E^ = I and 0 < E^, Eg < I, we 
can define a relative degree of excitation of particles in 
the two layers, by the threshold value 0.5: large exci¬ 
tations correspond to E^ g > 0.5, small excitations to 
E\ g < 0.5. If E\ = E’^ = 0.5, particle vibrations in 
both layers are of the same degree and correlated. 

At small w = 2, A we find, particularly in the low-w 
region, a large fraction of vibrational states with E\ ~ 
Eg ~ 0.5. As w increases, in the high frequency re¬ 
gion uj > ujg^^, where ujg^^ ~ 13 is the high-frequency 
boundary in gsiuj), only particles in the light layer A 
vibrate {E\ ~ 1), whereas those in the heavy layer B 
are almost immobile, as indicated by Eg ~ 0. In this w- 
region, phonon propagation in the cross-plane direction 


is therefore almost completely suppressed. On the other 
hand, for w < ujg^^, particles pertaining to the heavy 
layer B show large vibrational amplitudes {Eg > 0.5), 
while vibrations in layer A tend to be small {E^ < 0.5). 
More in details, for w > 20, we see that the averaged 
amplitudes are much larger in the B layer {{Eg') > 0.8) 
than in the A layer ((A^) < 0.2) in the 2 < uj < 7.5 
range. Contrary to the case of a; > ujg’^^, however, a 
significant number of modes are excited in both layers A 
and B, even with E^ — Eg ~ 0.5. We therefore conclude 
that, for UJ < ujg^^, some phonons still propagate in the 
cross-plane direction, contributing to kcp- 

We note that our observation of the vibrational sepa¬ 
ration in both the vDOS and vibrational amplitudes is 
consistent with results reported previously [28, 43, 53]. 
Indeed, the simulation work of Ref. [28] reported a sepa¬ 
ration in the vDOS of the Si isotopic-superlattice (^®Si- 
"‘^Si superlattice). A recent simulation work [43] focused 
on partial inverse participation ratios in a superlattice 
similar to the one considered here, reporting vibrational 
modes separation between layers. Ref. [53] attributed 
the reduction of thermal conductivity to a mechanism 
described as phonon localization, which we consider to 
be essentially the same phenomenon as the vibrational 
separation described here. 

We believe that this concept of vibrational separation 
is a simple and accurate framework to rationalize the 
behaviour of thermal conductivity in superlattices. In 
particular, it provides a complete characterization of 
the minimum in the w-dependence of kcp- Indeed, 
in the range w = 2 to 20 identifying the coherent 
regime, the vibrational separation hinders the coherent 
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FIG. 5. Thermal conductivity in superlattice S2 com¬ 
posed of ordered crystalline layers intercalated with 
mass disordered layers. The two components kcp and kip 
are plotted as functions of w. The mass ratio of the disor¬ 
dered alloy layer mB 2 /mBi is (a) 2, (b) 4, and (c) 8. The 
values Kcp kJp of the diffuse limits, together with those 
in the glass and the disordered alloy are indicated by lines. 


phonon propagation in the cross-plane direction, leading 
to the large reduction of kcp- In contrast, in-plane 
phonon propagation is very mildly affected by the 
vibrational separation and, therefore, Kip keeps high 
values. Also, by considering (A^) and (Eq) (solid 
lines), we conclude that the separation saturates to 
its maximum level at w ~ 20. Upon further increase 
w > 20, although averaged values show no significant 
changes, we recognize an increasing fraction of modes 
with > 0.5 and Eg < 0.5 for uj < (panels (e) 
w = AQ and (f) ic = 80 in Fig. 4). This observation 
indicates that the separation tendency for modes with 
E\ < 0.5 and Eg > 0.5 becomes weaker, i.e., the 
correlation of vibrational features in the two layers 
decreases, which corresponds exactly to the incoherent 
transport picture, and leads to the increase of kcp- 
Although transport becomes completely incoherent only 
for values of w of the order of the Kapitza length (note 
that £z ~ 1600 for rriBlmA = 4), this feature appears as 


soon as the vibrational separation is saturated, at the 
crossover point w* ~ 20. Thus, the saturation point of 
the vibrational separation identifies the minimum value 
of KCP) which can be indeed below the glass limit. 

S2. Superlattice composed of intercalated or¬ 
dered crystalline layers and mass disordered lay¬ 
ers. This system consists of three components, with 
masses = 1 in the crystalline layer A, and msi 
and mB 2 in the disordered alloy layer B. In Fig. 5, 
we plot KCP and kip for the mass ratios of the layer B, 
mB 2 /mBi = 2, 4, and 8. At small w < A, the values of 
both KCP and Kip are very close to those of the disordered 
bulk alloy formed by the same particles. As w increases. 
Kip increases gradually toward . This increase is con¬ 
trolled by the development of in-plane phonon propaga¬ 
tion in the ordered crystalline layer A. Indeed, the g{uj) 
of the superlattice, shown in Fig. 6, roughly follows that 
of the disordered bulk alloy at small w = A, whereas at 
large w = 20,40 it is dominated by gA{uj). In particular, 
the longitudinal peak around w ~ 14.5 becomes clear. 



0 5 10 15 20 25 


UJ 

FIG. 6. Vibrational density of states in superlattice 

S2. We show our results for the mass ratio of the disordered 
alloy layer mB^lmBi = 4, with msi = 0.4 and mB 2 = 1.6. 
The period w is 4, 20, and 40 for (a), (b), and (c), respectively. 
For comparison, we plot gA{u}) of the bulk crystal formed by 
particles of mass rriA ~ 1 (layer A), gaiuj) of the disordered 
alloy with masses msi = 0.4 and mB 2 = 1.6 (layer B), and 
the vDOS of the glass and the disordered alloy formed by the 
same constituent species. 
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W (monolayers) 

FIG. 7. Thermal conductivity in superlattice S3 com¬ 
posed of identical crystalline layers with weakened in¬ 
terface. Thermal conductivities kcp and kip are plotted as 
functions of w. The interface interaction €ab is 0.5 in (a) and 
0.1 in (b). We also show, by the horizontal lines, the thermal 
conductivities of the corresponding one-component homoge- 
neons bnlk crystal and glass with nnmodified interactions. 


corresponding to that of the crystalline layer A. 

The cross-plane value kcp also increases with w, but 
reaches the limit value n'^p already at w ~ 20. Since 
Kb of the disordered alloy layer B is low (see Table I), 
K^P remains low, typically less than twice the disordered 
alloy value. As a result, the variation of kcp with w is 
small. This result indicates that scattering in the disor¬ 
dered alloy layer B dominates the thermal conduction in 
the cross-plane direction. Both experimental work [54] 
on Si(crystal)-SiGe(disordered alloy) nanowires and nu¬ 
merical simulations [42] have reported similar observa¬ 
tions. We also note that the coherent nature of the su¬ 
perlattice phonons in the cross-plane direction, which we 
observed in the SI system, breaks down in S2. This is es¬ 
sentially equivalent to previous findings that disorder in 
interfacial roughness [39-41], or interfacial species mix¬ 
ing [42, 43] destroy the coherent features of vibrational 
excitations present in the investigated superlattices. In 
addition, the thermal conductivity tensor becomes in¬ 
creasingly anisotropic at larger w due to the increase of 
Kip , showing a behaviour different than that observed in 
SI. As a consequence of these features, in superlattices 
of type S2 the variability of the cross-plane heat trans¬ 
port is strongly bounded, and the minimum limit of kcp 
just corresponds to the disordered alloy limit, i.e., kcp 
cannot be reduced below the glass limit. 

S3. Superlattice composed by identical crys¬ 


talline layers separated by weakly interacting in¬ 
terfaces. In Fig. 7 we show the w-dependences of kcp 
and Kip for the case where the energy scale associated 
to particles interactions across the interfaces (sab) are 
lowered compared to those intra-layers, with eab = 0.5 
and 0.1 in the two panels. In the figure, we also plot as 
lines the data for the corresponding one-component crys¬ 
tal and glass with unmodified interactions. The in-plane 
value Kip is almost independent of w, and is very close 
to the value pertaining to the crystal. In contrast, kcp 
decreases monotonically by increasing w, and especially 
in the weaker case cab =0 .1, the observed reduction of 
KCP is dramatic. At w = 10, kcp equals the value ob¬ 
tained for the glassy sample, and it is almost two orders 
of magnitude lower than this value at w = 20. This ex¬ 
tremely low KCP is consistent with previous experimental 
work [21]. 

Some insight about the origin of this observation comes 
from the data shown in Fig. 8, where we display the aver¬ 
age cross-plane distance 5z between adjacent crystalline 
planes (monolayers), normalized to the value in the per¬ 
fect lattice, a/2. For cab = 0.5 and w = i, the system 
keeps the perfect lattice structure, with Sz = a/2 for all 
monolayers. In contrast, as cab decreases and for a large 
value w = 20, Sz becomes substantially larger than a/2 
at the interfaces, which therefore assumes a local density 
lower than the average. At the same time, slightly re¬ 
duced Sz are also observed for the other intra-monolayers, 
leading to an increase of the local density compared to the 
average. This heterogeneity hinders energy propagation 
across the interface and, as a result, phonons are spec¬ 
ularly reflected and confined in the in-plane direction. 
We remark that in the cases with w = 20, the values of 
Sz at the interfaces located at w/2 = 10 and u> = 20 
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FIG. 8. Distance between adjacent monolayers in su¬ 
perlattice S3. The average cross-plane distance Sz between 
adjacent crystalline planes plotted for each monolayer, iden¬ 
tified by the corresponding order index. We present the value 
of Sz normalized to a/2, the horizontal line Sz/{a/2) = 1 
therefore indicates the value in the perfect crystalline lattice. 
The displacements observed in the cases w = 20 are discussed 
in the main text. 
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FIG. 9. Vibrational density of states and vibrational amplitudes in superlattice S3. We report data corresponding 
to the indicated values of the the interfacial interaction energy cab and the repetition period w. (a) cab = 0.5, w = 4, (b) 
CAB = 0.5, w = 20, and (c) cab = 0.1, w = 20. In panels at the top, we show the vDOS for the superlattices of type S3, 
together those of the corresponding one-component homogeneous crystal and glass with unmodified interactions. In panels at 
the bottom we show the vibrational amplitudes, E'X and E%, in layers A and B, plotted as functions of the eigenfrequency oj*. 
The solid line represents the average values (E'X) and {E%), calculated in bins of the form ± Soj'^/2, with 6uj'° = 0.5. The 
horizontal dotted lines indicate i?A = = 0.5. 



are different, with a large discrepancy for eab = 0.1. 
We rationalize this behaviour by observing that, during 
the preparation stage of the sample, the applied selec¬ 
tive weakening of the interactions destabilizes the global 
equilibrium of the superlattice, with a concentration of 
mechanical stress close to the interfaces. Lattice planes 
far from the boundaries easily recover mechanical equi¬ 
librium by coherently reducing their mutual distance. In 
contrast, particles in monolayers adjacent to the inter¬ 
faces move both out-of-plane and in-plane, to optimize 
the local effective spring constants. The optimal solution 
found depends in general on the details of the local en¬ 
vironment, explaining the observed discrepancy in Sz at 
different interfaces. 

The behaviour of kcp can be further elucidated by in¬ 
spection of the main features of the vibrational spectrum. 
In Fig. 9 we plot the g{uj) of superlattice S3, together with 
the vibrational amplitudes E\ and Eg. The g{uj) shows 
transverse and longitudinal phonon branches for all cases, 
similar to the homogeneous bulk crystal. As w increases 
g{u}) deforms, following the appearance of an increasing 
fraction of modes at increasing higher frequencies. This 
behaviour is certainly correlated to the observation made 
above (see Fig. 8) for w = 20, that the distance between 
monolayers far from the interfaces becomes smaller than 
a/2. The consequent larger mass density makes higher 
the frequency of phonon modes of given wavelength, lead¬ 
ing to the shift of g{uj) towards higher frequencies. This 
global shift has as a consequence a mild increase of kip 
with w, as it is clear from Fig. 7(b) (cab = 0.1 case). 
Note that for cab = 0.5 and w = 4 (Fig. 9(a)), g{uj) 
shows an excess of lower-w modes compared to those 


present in the one-component crystal, simply due to the 
weakened interactions at the interfaces. 

We now focus on the vibrational amplitudes, E\ and 
Eg (Fig. 9, bottom panels). In the cases with cab = 0.5 
and u> = 4 and 20, the particles in the two layers A 
and B show completely equivalent and correlated vibra¬ 
tions for the vast majority of the modes, as indicated 
by E\ = Eg = 0.5. This result implies that phonons 
indeed propagate across the weakened interfaces in the 
cross-plane direction, but they are also partially reflected 
at the interface, causing the observed reduction of kcp- 
The situation changes drastically in the case cab = 0.1 
and w = 20, where the ultra-low value of kcp can be 
reached. Except for the low-w modes, E\ and Eg are 
symmetrically randomly distributed around the average 

values = 0.5, indicating that particles in layers 

A and B vibrate independently, in an uncorrelated man¬ 
ner. As a consequence, a very large fraction of vibrational 
modes do not cross at all the interfaces, but rather un¬ 
dergo a perfect specular reflection. In this situation, heat 
is not transferred between two adjacent layers A and B, 
leading to extremely low value of kcp, while keeping a 
high Kip. We conclude by noticing that although spec¬ 
ular reflection was also observed in the system SI, the 
physical mechanism behind this phenomenon is different 
in the two cases: vibrational separation causes reflection 
in the former, whereas weakened interactions across the 
interface, with the resulting augmented spacing between 
the layers, completely block cross-plane phonon propa¬ 
gation in the latter. 
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III. DISCUSSION 


We have provided numerically, for the first time to 
our knowledge, a clear demonstration of very low ther¬ 
mal conductivities in super lattices, below the glassy limit 
of the corresponding amorphous structures. Blocking 
phonon propagation in ordered structures via interfaces 
design is the key principle. We have identified two possi¬ 
ble strategies to achieve this goal: imposing a large mass 
heterogeneity in the intercalated layers (as in system SI) 
or degrading inter-layers interactions compared to those 
intra-layers (as in S3). We have found that in both cases 
phonons are specularly reflected at the interface and con¬ 
fined in the in-plane direction. This reduces the cross¬ 
plane thermal conductivity kqp below the corresponding 
glass limit, while keeping the in-plane contribution kip 
close to the pure crystalline value. 

More specifically, in the case of mass mismatch (SI), 
propagation of phonons with high frequencies (w > 
is almost completely suppressed, whereas a fraction of 
low-frequency phonons (w < are still able to 

propagate across the interfaces, contributing to kcp 
(F ig. 4(d)). Also, the minimum in thermal conductivity 
as a function of the repetition period w (Fig. 2) cor¬ 
responds to a maximum in the vibrational separation 
between the layers of type A and B. These therefore 
act as true filters in complementary regions of the vibra¬ 
tional spectrum, suppressing significantly phonons trans¬ 
port in the direction of the replication pattern. On the 
other hand, attenuated interactions across the interfaces 
{S3) are able to block phonons at almost all frequencies 
(see Fig. 9(c)), which results into extremely low values 
of KcP) even orders of magnitude lower than the corre¬ 
sponding glass limit (Fig. 7(b)). In this sense, directly 
modifying the interfaces seems to be the most effective 
strategy to obtain very low heat transfer. Note that 
this is a practically feasible route, since attenuated in¬ 
terfaces can be designed by exploiting materials with 
weak van der Waals forces among adjacent crystalline 
planes, as demonstrated in the case of WSe 2 sheets in 
Ref. [21]. Interfaces stiffness modification by controlling 
pressure [46, 47] or chemical bonding [48] are additional 
possible routes to directly tune the strength of interfaces. 

Our data also suggest that intercalating disordered al¬ 
loy layers in ordered crystalline layers {S2) is not effective 
in lowering kcp- Indeed, we have demonstrated that in 
this case disorder is not sufficient to block the propa¬ 
gation of vibrational excitations, even though it makes 
phonons lifetimes short. The intercalated disordered al¬ 
loy layer dominates phonon transport in the entire super¬ 
lattice, notwithstanding the presence of the crystalline 
layers. As a result, thermal conductivity is very similar 
to the one of the disordered alloy and is only marginally 
modified by modulation of the period w (see Fig. 5). 
Also, as suggested in previous works, disorder in the in¬ 
terfacial roughness [39-41] or interfacial mixing [42, 43] 
seems to already dominate over phonon transport, and 
destroy the coherent nature of phonons. 


In addition, as we understand from our analysis of vi¬ 
brational amplitudes (Figs. 4 and 9), it is much more 
problematic to block low-w (long wavelength. A) phonons 
propagation, than those with high-w (short A). This situ¬ 
ation is similar to what has been observed in bulk glasses, 
where the long-A acoustic waves are not scattered by the 
disorder and can propagate over long distances by car¬ 
rying heat energy [55, 56]. Therefore, blocking or effi¬ 
ciently scattering the long-A phonons is also a key fac¬ 
tor to achieve very low thermal conductivities, as was 
pointed out in Ref. [57]. A possibility to realize this task 
is embedding in the targeted material objects featuring 
larger typical sizes, including nano-particles [58, 59] or 
nano(quantum)-dots [22, 60]. Based on this strategy, 
very low thermal conductivity was achieved experimen¬ 
tally in a Si-Ge quantum-dot superlattice [22], even be¬ 
low the amorphous Si value. The additional possibility of 
introducing large size defects by the porous structuring 
of materials has also been explored in a recent numerical 
work [61]. Here, values of thermal conductivity 10"^ times 
smaller than that of bulk Si were reached in Si phononic 
crystals with spherical pores. 

In conclusion, we note that the three superlattice struc¬ 
tures studied in the present work show totally different 
zc-dependences of cross and in-plane thermal conductivi¬ 
ties. Our results therefore not only contribute to a deeper 
comprehension of the physical mechanisms behind very- 
low thermal conductivity, they also provide insight for 
developing new design concepts for materials with con¬ 
trolled heat conduction behaviour. 


IV. METHODS 


System description. In this Section we provide de¬ 
tails on the numerical models we have used for the su¬ 
perlattices. The corresponding amorphous structures 
(glasses) and disordered alloys with exactly the same 
composition were also prepared, for the sake of compar¬ 
ison with superlattice phases. We have considered in all 
cases a 3-dimensional cubic box, of volume V = {L 
being the linear box size), with periodic boundary condi¬ 
tions in all directions. In the superlattice and disordered 
alloy cases, particles were distributed on the FCC lattice 
sites. In the glass phases, they were frozen in topolog¬ 
ically random positions following a rapid quench from 
the normal liquid phase below the glass transition tem¬ 
perature Tg, avoiding crystallization (see, for instance. 
Ref. [55] for details on the preparation of glasses). Par¬ 
ticles, i and j, interact via soft-sphere (SS) or Lennard- 
Jones (LJ) potentials: 


v^s(r)=e^^ 







( 3 ) 


where r is the distance between those two particles, and 
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a’’^ and e*-' are the interparticle diameter and interaction 
energy scale, respectively. The potential is cut-off and 
shifted at Tc = 2.her*-’. Particle i has mass m*, and we 
have used a, e/ks (ks is the Boltzmann constant), and m 
as units of length, temperature, and mass. As a reference, 
for Argon a = 3.4A, e/fcs = 120 K, and m = 39.96 a.m.u. 
We considered the number density p = N/V = 1.015, 
corresponding to a lattice constant a = (4/p)^/^ = 1.58. 

We prepared three superlattices, composed of inter¬ 
calated FCC lattice layers, A and B, both of thickness 
w/2, as schematically illustrated in Fig. 1. The first su¬ 
perlattice (SI) consists of two crystalline layers formed 
by sphere particles with different masses, uia and niB- 
We have considered mass ratios uib/ttia > 1, while keep¬ 
ing a constant average mass {itia + mB)/^ = 1. As an 
example, the case niB/mA = 4 corresponds to uia = 0.4 
and niB = 1.6. We have dubbed A and B as the light 
and heavy layers, respectively. Note that a mass ratio 
of TTiBlmA = 2.5 corresponds to the case of the realistic 
Si-Ge superlattice. Except for the above mass difference 
in the different layers, all particles are characterized by 
the same properties. In particular, they interact via the 
SS potential Ugg(r), with = 1. 

The second superlattice {S2) is composed of an ordered 
crystalline layer A intercalated to a disordered alloy layer 
B. TTiA = 1 in A, whereas in B half of the particles have 
mass niBi, 'rnB2 the others, and are randomly distributed 
on the lattice sites. Again, tubi and mB 2 are determined 
by the mass ratio mB 2 /iTiBi > Ij keeping a constant 
average value (itibi + rnB2)/‘^ = 1- All particles in both 
layers interact via the SS potential Vgg(r) with ct*-! = 
= 1 . 

The third superlattice {S3) is composed of identical 
crystalline layers A and B, but the interactions among 
particles in different layers (i.e., across the interfaces) 
are modified (weakened) compared to those intra-layers. 
All particles have mass niA = ttib = 1, and interact 
via the LJ potential v^ij{r), with cr*-! = e*-! = 1. The 
energy scale of interactions between particles pertaining 
to different layers are, however, reduced to = cab < 1 - 

MD simulation and the Green-Kubo method for 
the calculation of thermal conductivity. In the 

present study, all simulations have been realized by us¬ 
ing the large-scale, massively parallel molecular dynam¬ 
ics simulation tool LAMMPS [62, 63]. The systems were 
first equilibrated at relatively low temperature T = 10“^ 
by MD simulation in the NVT-ensemh\e. This choice 
was dictated by the need to reduce anharmonic effects, 
in order to primarily focus on the contribution of the 
structural features of the superlattices on thermal con¬ 
ductivity. We must note that our approach is classical, 
and does not take into account the quantum mechanisms 
active in the low-T regime [10]. These effects have im¬ 
portant implications, increasing the contribution to the 
thermal conductivity coming from low-w vibrational ex¬ 
citations. At present, however, it is not obvious and still 
under debate how to effectively include quantum effects 


into a classical system [64, 65], and we have therefore 
chosen to stay within a fully classical approach. 

Following the equilibration stage, we performed the 
production runs in the NVE-ensemble. The Green-Kubo 
formulation [50, 51] was next applied to calculate the 
thermal conductivities, in the cross-plane and in-plane 
directions, respectively: 


KCP 

Kip 


^ {J,{t)JM)dt, 

2yj'2 J {Jx{t)Jx{b) + Jy{t)Jy{0)) dt. 


(4) 


Here, Jx,yi and Jz are the heat currents in the in-plane 
{x, y) and cross-plane (z) directions, and () denotes the 
ensemble average. In the bulk glasses and disordered 
alloys, Kcp — Kip, i.e., heat conduction is isotropic, 
whereas in the superlattices, they are expected to assume 
different values [18, 19]. Landry et al. [51] have carefully 
confirmed the validity of the Green-Kubo method for the 
calculation of superlattices thermal conductivity, by com¬ 
parison with the direct method based on non-equilibrium 
simulation. Also, in the Green-Kubo calculations, one 
must be attentive to finite system size effects [50, 51]. In¬ 
deed, long-wavelength phonons with X > L are excluded 
from the simulation box due to the finite value L of the 
box size, which imposes important size effects on the nu¬ 
merical determination of k. The box size therefore needs 
to be large enough to include a vibrational spectrum suffi¬ 
cient to establish an accurate description of anharmonic 
coupling (scattering) processes [50]. We note that the 
considered T = 10“^ is low enough to substantially re¬ 
duce anharmonic effects, but anharmonic couplings are 
still present. 

We can take care of finite size effects by increasing L 
to values where kcp and Kip become L-independent. For 
the glass and disordered alloy thermal conductivities, we 
have confirmed that a system size L = 10a {N = 4,000) 
is sufficiently large to obtain correct values of kqp — Kip, 
without any size effect [8, 9]. In the super lattice cases, 
the appropriate L depends on the considered structure 
and the periodic repetition length w [51]. More in de¬ 
tails, we paid particular attention to the number P of 
repetitions, defined from L = Pw, necessary to pro¬ 
duce sufficient anharmonic couplings of phonons in the 
cross-plane direction. We have therefore investigated the 
presence of finite size effects by analyzing different sys¬ 
tems with sizes ranging from L = 10a (20 monolayers, 
N = 4,000) to 24a (48 monolayers, N = 55,296). In 
Figs. 2, 5, and 7, we show multiple data points at some 
w-values, obtained for different system sizes. For the SI 
super lattice, we confirmed that the required number P of 
repetitions becomes larger for smaller w [51]: one period 
{L = w) only is adequate for w > 20, whereas four peri¬ 
ods or more {L > Aw) are required for w < 8. We have 
therefore employed four pattern repetitions {L = Aw) for 
10 < lu < 12 and two {L = 2w) for 14 < w < 18. 
This behaviour is simple to rationalize by inspecting the 
data in Fig. 2, where the crossover between incoherent 
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and coherent phonon transport occurs around w* ~ 20. 
In the coherent regime w < 20, the wave character of 
the phonons becomes important, and therefore a larger 
number of repetitions is necessary to produce the coher¬ 
ent wave interference processes correctly. In contrast, 
smaller values of P are needed (even P = 1) in the in¬ 
coherent regime w > 20, where the incoherent particle 
nature of the phonons appears. 


For the S2 and S3 superlattices the system size 
effects issue is much less pronounced than in the SI 
case. We can understand this behaviour by noticing 
that phonon tranport is mainly determined by the 
scattering processes in the disordered alloy layer in 
S2, and the blocking at the weak interface for S3. In 
both cases the missing long wavelength phonons, with 
\ > L, play very little role in phonon transport and 
finite system size effects are consequently negligible. We 
therefore used P = 1 {L = w) for w > 20 and one or 
more repetitions (L > w) for w < 20, for both S'^and S3. 


Normal modes analysis. We have characterized 
the superlattice vibrational states (superlattice phonons) 
by performing a standard normal-mode analysis [52] 
with ARPACK [66]. We have diagonalized the dy¬ 
namical (Hessian) matrix calculated at local minima 
of the potential energy landscape, and obtained eigen¬ 
values A* and eigenvectors (polarization vectors) = 
{ej,..., e^,..., e^}. Here, j is the particle index, and 
k = 1,2,..., SN — 3 is the eigenmode number, where we 
have disregarded the three vanishing Goldstone modes. 
The eigenvectors are normalized such that • e* = 
= ^ki, where 6ki is the Kronecker delta 
function. The eigenfrequencies are next calculated as 
u)^ = and the associated probability distribution 


(normalized histogram) directly provides the vDOS: 


5(w) 


I 

3N- 3 


3N-3 


( 5 ) 


In addition, from the eigenvector we have defined 
the vibrational amplitudes of mode k for layers A and B: 

Emb)= E ( 6 ) 

jelayer.4(B) 


Note that -I- = I for each k and, 

therefore, 0 < E\,Eg < I. Based on the values of 
and Eg, one can determine in which layer particles are 
more displaced (excited) according to the associated 
eigenvector e'^. In particular, if E\ > 0.5, Eg < 0.5 
{E\ < 0.5, Eg > 0.5), particles in layer A {B) contribute 
more to mode k than those in layer B (A). In the case 
E\ = Eg = 0.5, particles in both layers contribute 
equivalently, and in a correlated manner. Note that 
the normal mode analysis provides us with the system 
vibrational states in the harmonic limit T —>■ 0 which, 
we believe, is an appropriate approximation for our case 
T = 10“^, where anharmonicities are weak. 
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